home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Turnbull China Bikeride
/
Turnbull China Bikeride - Disc 1.iso
/
ARGONET
/
PD
/
MATHS
/
RLAB
/
RLAB125.ZIP
/
!RLaB
/
toolbox
/
gamma
< prev
next >
Wrap
Text File
|
1994-02-21
|
5KB
|
218 lines
//-------------------------------------------------------------------//
//
// Syntax: gamma ( A )
// gammaln ( A )
// gammp ( A , X )
// gammq ( A , X )
// Description:
// gamma: Computes the value of the gamma function. The input, A
// must be a row, or column vector. Allowable values of A
// are: A > 0.
// gammln: Computes the value of the log-gamma function.
// gammp: Computes the value of the incomplete gamma function,
// P(a,x).
// Reference: Numerical Recipes in C, pp 172.
// gammq: Computes the value of the incomplete gamma function,
// Q(a,x) = 1 - P(a,x).
// Reference: Numerical Recipes in C, pp 173.
//-------------------------------------------------------------------//
static (gser, gcf);
gamma = function( A )
{
if (!all (all (A >= 0))) { error ("valid range for gamma() 0:inf"); }
return exp( gammaln( A ) );
};
gammaln = function( A )
{
local(a, cof, j, ser, tmp);
if( !all (all (A >= 0)) ) // Test for correct input range
{
return inf();
}
cof = [76.18009173, -86.50532033, 24.01409822, ...
-1.231739516, 0.120858003e-2, -0.536382e-5];
a = A - 1;
tmp = a + 5.5;
tmp = tmp - (a + 0.5) .* log(tmp);
ser = 1;
for( j in 1:6 )
{
a = a + 1;
ser = ser + cof[j]./a;
}
return ( -tmp + log(2.50662827465*ser) );
};
//
// Computes the value of the incomplete gamma function, P(a,x).
// Currently the input, A, must be a row, or column matrix.
//
gammp = function( X, A )
{
local(a, gammcf, gammser, gln, g, i);
if( any (any (X < 0)) || any (any (A < 0))) {
error ("invalid arguments to gammp()");
}
if (!all (size (X) == size (A)))
{
if (A.n != 1)
{
error ("size of 2nd arg must match 1st, or be 1-by-1");
else
a = A.*ones (size (X));
}
else
a = A;
}
g = zeros (size (X));
for (i in 1:X.n)
{
if(X[i] < (a[i] + 1))
{
g[i] = gser (a[i], X[i], gln);
else
g[i] = 1 - gcf (a[i], X[i], gln);
}
}
return g;
};
//
// Computes the value of the incomplete gamma function,
// Q(a,x) = 1 - P(a,x).
//
gammq = function( X , A )
{
local(a, gammcf, gammser, gln, g, i);
if(any (any (X < 0)) || any (any (A <= 0))) {
error ("invalid arguments to gammq()");
}
if (!all (size (A) == size (X)))
{
if (A.n != 1)
{
error ("size of 2nd arg must match 1st, or be 1-by-1");
else
a = A.*ones (size (X));
}
else
a = A;
}
g = zeros (size (X));
for (i in 1:X.n)
{
if(X[i] < (a[i] + 1))
{
g[i] = 1 - gser ( a[i], X[i], gln );
else
g[i] = gcf ( a[i], X[i], gln );
}
}
return g;
};
//-------------------------------------------------------------------//
//
// Syntax: gser ( A, X, GLN )
// Description:
// Computes the value of the incomplete gamma function, evaluated by
// it's series representation.
// Reference: Numerical Recipes in C, pp 173.
//-------------------------------------------------------------------//
gser = function( A, X, GLN )
{
local(ITMAX, ap, del, eps, n, sum);
GLN = gammaln(A);
if( any (any(X <= 0)))
{
if( any( X < 0 ) ) {
error("X less than 0 in function gser()");
}
return 0;
else
ITMAX = 100; eps = epsilon();
ap = A;
del = sum = 1./A;
for( n in 1:ITMAX )
{
ap = ap + 1;
del = del .* X./ap;
sum = sum + del;
if( all (all(abs (del) < abs (sum).*eps)))
{
return sum.*exp(-X + A.*log(X) - GLN);
}
}
error("A too large, ITMAX too small in function gser()");
return 0;
}
};
//-------------------------------------------------------------------//
//
// Syntax: gcf ( A, X, GLN )
// Description:
// Computes the value of the incomplete gamma function, evaluated by
// it's continued fraction representation.
// Reference: Numerical Recipes in C, pp 174.
//-------------------------------------------------------------------//
gcf = function( A, X, GLN )
{
local(ITMAX, a0, a1, an, ana, anf, b0, b1, eps, fac, g, gold, n);
ITMAX = 100;
a0 = 1; b0 = 0; b1 = 1; fac = 1; gold = 0;
eps = epsilon();
GLN = gammaln(A);
a1 = X;
for( n in 1:ITMAX )
{
an = n;
ana = an - A;
a0 = (a1 + a0.*ana).*fac;
b0 = (b1 + b0.*ana).*fac;
anf = an.*fac;
a1 = X.*a0 + anf.*a1;
b1 = X.*b0 + anf.*b1;
if( all(a1) )
{
fac = 1./a1;
g = b1.*fac;
if( abs((g - gold)/g) < eps )
{
return exp(-X + A.*log(X) - GLN).*g;
}
gold = g;
}
}
error("A too large, ITMAX too small in function gcf()");
};